close all
clear all
clc

Eng=13.474;
lambda=0.9201;
Attl=4854.22e4; % Attenuation length in unit of Angstrom for water @ 13.474 keV
mu=1/Attl;
beta=mu*lambda/4/pi;

I0=1e4;  

Fluores_Fe_Qz_file={
    'sph01_scan329_Fe_EL.txt';     % 1:  bar 3  SPM/KCl/FeCl3 pH 4
    'sph02_scan406_Fe_EL.txt';     % 2:  bar 3 KCl 100 mM, FeCl3 5 mM
    }

qz=[0.01:0.001:0.035]';
for k=1:1
    data=dlmread(Fluores_Fe_Qz_file{k});
    FeEL{k}.Intensity=data(:,2);
    FeEL{k}.Intensity_err=data(:,3);
end

data=dlmread(Fluores_Fe_Qz_file{2});
Intensity_sub=data(:,2);
Intensity_sub_err=data(:,3);


MarkerType={'o','s','v','^','d','>','<'}
COLOR={'k','r','b',[0 0.4 0],'m','k','k'}

%% Subphase only profile fitting
close all
clc
Param=[
    1, 0.0219, 0
    ]
LB=[0, 0.0218-1e-6,  0];
UB=[1, 0.022+1e-6,  1 ];

xfit=(0.01:0.0001:0.0300)';
xdata_sub=qz(Intensity_sub>0 & qz<0.031);
ydata_sub=I0*Intensity_sub(Intensity_sub>0 & qz<0.031);
ydata_sub_err=I0*Intensity_sub_err(Intensity_sub>0 & qz<0.031);
[p_opt_sub]=optim_fitsubphase(Param,xdata_sub,ydata_sub,ydata_sub_err,LB,UB)
C_bulk=p_opt_sub(1);
Qc=p_opt_sub(2);
bkg=p_opt_sub(3);
T=Transmission(xfit,Qc,beta,lambda);
D=fitsubphase(xfit,Qc,Eng, Attl);
yfit_sub=C_bulk*T.*D+bkg;
figure
hold on
H=ploterr(xdata_sub,ydata_sub,[],ydata_sub_err,MarkerType{1},'hhy',0.1);
set(H,'color',COLOR{1},'LineWidth',3,'MarkerSize',8);
plot(xfit,yfit_sub,'-','LineWidth',3,'color',COLOR{k});
hold off
set(gca,'Linewidth',3,'box','on');
set(gca,'FontSize',30,'FontName','times');
set(gca,'xlim',[0.009,0.036]);
% set(gca,'ylim',[-0.1 3.1]);
axis square

%%
close all
clc


figure
hold on
for k=1:1
    FeEL_corr{k}.Intensity=I0*(FeEL{k}.Intensity);
    FeEL_corr{k}.Intensity_err=I0*sqrt(FeEL{k}.Intensity_err.^2);
    H=ploterr(qz,FeEL_corr{k}.Intensity,[],FeEL_corr{k}.Intensity_err,MarkerType{k},'hhy',0.1);
    set(H,'LineWidth',3,'color',COLOR{k})
end
%%
close all
clc

figure
hold on
xfit=(0.01:0.0001:0.0300)';
qmin=1;
qmax=21;
C_EL=zeros(1,5);
Zion=zeros(1,5);
for k=1:1
    x=qz;
    y=FeEL_corr{k}.Intensity;
    yerr=FeEL_corr{k}.Intensity_err;
    
    xdata=x(qmin:qmax);
    ydata=y(qmin:qmax);
    ydata_err=yerr(qmin:qmax);
    
    
    %%%%%%%%%%%%%%%%%%%%% Radiation Damage Model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Param=[
%         1, 0.0218, 0, 0.1,0,0.021
%         ]
%     LB=[0, 0.0218,  20, 0, 0.022];
%     UB=[100, 0.0220, 100,1,0.035];
%     [p_opt{k}]=optim_fitsurface_tau(Param,xdata,ydata,ydata_err,LB,UB) 
%     C_EL=p_opt{k}(1);
%     Zion=p_opt{k}(3);
%     tau=p_opt{k}(4);
%     qtau=p_opt{k}(5);
%     
%     T=Transmission(xfit,p_opt{k}(2),beta,lambda);
%     D=fitsubphase(xfit,p_opt{k}(2),Eng, Attl);    
%     yfit{k}=C_EL*T.*exp(-Zion./D).*exp(-0.5*(1+erf((xfit-qtau)/0.0001))*tau.*(xfit-qtau)./qtau);
    
    %%%%%%%%%%%%%%%%%%% Regular Model for esitimation maximum binding %%%%%%%%%%%%%%%%%%%%
    Param=[
        1, 0.0218, 10
        ]
    LB=[0, 0.0218-1e-6,  10];
    UB=[100, 0.022+1e-6, 50];
    
    p_opt=optim_fitsurface(Param,xdata,ydata,ydata_err,LB,UB);
    C_EL(k)=p_opt(1);
    Qc=p_opt(2);
    Zion(k)=p_opt(3);
    T=Transmission(xfit,p_opt(2),beta,lambda);
    D=fitsubphase(xfit,p_opt(2),Eng, Attl);    
    yfit{k}=C_EL(k)*T.*exp(-Zion(k)./D);
    H=ploterr(x,y,[],yerr,MarkerType{k},'hhy',0.1);
    set(H,'color',COLOR{k},'LineWidth',3,'MarkerSize',8);
    plot(xfit,yfit{k},'-','LineWidth',3,'color',COLOR{k});
end
set(gca,'Linewidth',3,'box','on');
set(gca,'FontSize',30,'FontName','times');
set(gca,'xlim',[0.009,0.031]);
set(gca,'ylim',[-0.1 5.1]);
axis square


%% Calculate the surface density
% subphase concentratin is 5 mM equivalent to 5*
nb=5*1e-3*6.02e23/1e27;
NumFeDensity=zeros(5,1);
for k=1:5
    NumFeDensity(k)=C_EL(k)/C_bulk*nb;
end
NumFeDensity
Zion
concentration=[2, 5, 10, 20, 40];
close all
clc
figure
plot(concentration, NumFeDensity,'ko','MarkerSize',10,'LineWidth',3);
set(gca,'linewidth',3);
set(gca,'box','on');
set(gca,'FontSize',30);
set(gca,'xlim',[-1 50]);
axis square

% # of Molecule per surface area =  C_surf/C_bulk*nb
%% Display
close all
figure
hold on
plot(xfit,yfit_sub,'-','LineWidth',3,'color','k');
H=ploterr(xdata_sub,ydata_sub,[],ydata_sub_err,MarkerType{1},'hhy',0.1);
set(H,'color','k','LineWidth',3,'MarkerSize',14,'MarkerFaceColor','w');
plot(xfit,yfit{1},'-','LineWidth',3,'color','r');
H=ploterr(xdata,ydata,[],ydata_err,MarkerType{2},'hhy',0.2);
set(H,'color','r','LineWidth',3,'MarkerSize',14,'MarkerFaceColor','w');
hold off
set(gca,'xlim',[0.0095 0.0305],'ylim',[-0.5 5.5]);
set(gca,'fontsize',30,'linewidth',3);
set(gca,'box','on');
axis square


%% Error Analysis
% Comment: due to the potential radiation damage issue, the depth of the
% ion layer is considered to be the main error/uncertainty source
% So just shake up the Zion
close all
clc

Zion_test=[25:1:55];

UB_test=UB;
LB_test=LB;
Param=[
    1, 0.0218, 10
    ]
LB=[0, 0.0219-1e-6,  0];
UB=[100, 0.0219+1e-6, 50];

deltaZion=zeros(1,5);
ZionOpt=zeros(1,5);
chisq_min=zeros(1,5);
IronSurfDensity_opt=zeros(1,5);
IronSurfDensity_err=zeros(1,5);
 
for filenum=1:1
    x=qz;
    y=FeEL_corr{filenum}.Intensity;
    yerr=FeEL_corr{filenum}.Intensity_err;
    
    xdata=x(qmin:qmax);
    ydata=y(qmin:qmax);
    ydata_err=yerr(qmin:qmax);
    
    Z=[];
    chisq=[];
    IronSurfDensity=[];
    for k=1:length(Zion_test)
        UB_test(3)=Zion_test(k)+1e-6;
        LB_test(3)=Zion_test(k)-1e-6;
        [p_test,fval]=optim_fitsurface(Param,xdata,ydata,ydata_err,LB_test,UB_test);
        
        Z=[Z,Zion_test(k)];
        chisq=[chisq,fval];
        IronSurfDensity=[IronSurfDensity, p_test(1)/C_bulk*nb];
        
        T=Transmission(xfit,p_test(2),beta,lambda);
        D=fitsubphase(xfit,p_test(2),Eng, Attl);
        yfit{k}=p_test(1)*T.*exp(-Zion_test(k)./D);
      
    end
   
    figure
    subplot(1,3,1);
    hold on
    coef=polyfit(Z,chisq,2);
    chisq_min=coef(3)-coef(2)^2/coef(1)/4;
    chisqfit=polyval(coef,Z);
    plot(Z,chisq,'ko','MarkerSize',10);
    plot(Z,chisqfit,'k-','linewidth',3);
    hold off
    deltaZion(filenum)=sqrt(1/coef(1))*sqrt(chisq_min) % Rescale the uncertainty in parameter.
    ZionOpt(filenum)=-coef(2)/coef(1)/2
    set(gca,'box','on','LineWidth',3);
    grid on
    subplot(1,3,2);
    hold on
    plot(Z,IronSurfDensity,'bs','MarkerSize',10);
    %IronSurfaceDensity_opt(filenum)=IronSurfDensity( ZionOpt(filenum)-deltaZion(filenum)<Z<ZionOpt(filenum)+deltaZion(filenum))
    coef=polyfit(Z,IronSurfDensity,1);
    IronSurfDensityFit=polyval(coef,Z);
    plot(Z,IronSurfDensityFit,'k-','linewidth',3);
    IronSurfDensity_opt(filenum)=polyval(coef,ZionOpt(filenum))
    IronSurfDensity_err(filenum)=coef(1)*deltaZion(filenum)
    hold off
    set(gca,'box','on','LineWidth',3);
    grid on
   
    subplot(1,3,3);
    hold on
    H=ploterr(xdata,ydata,[],ydata_err,'ko');
    set(H,'MarkerSize',10,'LineWidth',3);
    for k=1:length(Zion_test)
        plot(xfit,yfit{k},'k-','LineWidth',3);    
    end
    hold off
    set(gca,'xlim',[0.009,0.036]);
    set(gca,'box','on','LineWidth',3);
    grid on
end

%% plot error analysis results
close all
clc
figure
subplot(1,2,1)
%hold on
H=ploterr(concentration, IronSurfDensity_opt,[],IronSurfDensity_err,'ko','hhy',0.1)
set(H,'MarkerSize',10,'LineWidth',3);
set(gca,'linewidth',3);
set(gca,'box','on');
set(gca,'FontSize',30);
set(gca,'xlim',[-1 50]);
axis square


% % %%
% % 
% % haxes1=gca;
% % 
% % haxes1_pos = get(haxes1,'Position'); % store position of first axes
% % haxes2 = axes('Position',haxes1_pos,...
% %               'XAxisLocation','top',...
% %               'YAxisLocation','right',...
% %               'Color','none');
% %           
% % H=ploterr(concentration,ZionOpt,[],deltaZion,'bs','hhy',0.1)
% % set(H,'MarkerSize',10,'LineWidth',3,'parent','haxes2');
% % % set(gca,'linewidth',3);
% % % set(gca,'box','on');
% % % set(gca,'FontSize',30);
% % % set(gca,'xlim',[-1 50]);
% % % set(gca,'ylim',[-1 51]);
% % % set(H,'yaxislocation','right');
% % axis square
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
% % 
